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Abstract: Inspired by recent results in the context of AdS/CFT integrability, we recon¬ 
sider the Thermodynamic Bethe Ansatz equations describing the ID fermionic Hubbard 
model at finite temperature. We prove that the infinite set of TBA equations are equivalent 
to a simple nonlinear Riemann-Hilbert problem for a finite number of unknown functions. 
The latter can be transformed into a set of three coupled nonlinear integral equations 
defined over a finite support, which can be easily solved numerically. We discuss the 
emergence of an exact Bethe Ansatz and the link between the TBA approach and the 
results by Jiittner, Kltimper and Suzuki based on the Quantum Transfer Matrix method. 
We also comment on the analytic continuation mechanism leading to excited states and 
on the mirror equations describing the finite-size Hubbard model with twisted boundary 
conditions. 
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1. Introduction 


The Hubbard model [1] arises as an approximate description of correlated electrons in 
narrow-band materials. Despite its simplicity, it is believed to capture important nonper- 
turbative features of real many-body fermionic systems. In one dimension, the model is 
exactly solvable and its relevance for the study of strongly correlated electrons is perhaps 
comparable to that of the Ising model for magnetism. Quite unexpectedly, the Hubbard 
model has also attracted the attention of high energy physicists due to its multiple con¬ 
nections [2, 3, 4, 5, 6, 7, 8] with the integrable spin chains emerging in the context of A7=4 
Super Yang-Mills theory [9, 10, 11]. In particular, the exact Bethe Ansatz equations of 
Lieb and Wu [12], and the corresponding string hypothesis [13], have many features in 
common with the asymptotic Bethe Ansatz of Beisert and Staudacher for the anomalous 
dimensions of single-trace operators with large quantum numbers [11]. 

We start from the ID Hubbard Hamiltonian written in the form 
L L 

H 0 = "E E ( C la c i+ 1,0- + c !+1,<t C v) + u - 2c l,t C ht)( 1 - ( L1 ) 

i=1 <7=t,t i =1 

where L is the length of the chain, c| and Cj. CT are fermionic creation-annihilation operators 
satisfying 

{ C i,cr> Cj,r} = { c J jC r> C j,ri = c y r } = (1-2) 

with periodic boundary conditions c\ A = cl+i,<t, c{ a = c* L+1 CT , and u > 0 is a dimensionless 
coupling constant proportional to the electric charge. For L even, the Hamiltonian (1.1) is 
50(4)-symmetric [14]. We shall consider a two-parameter deformation of (1.1) where this 
symmetry is explicitly broken by coupling the electrons to a chemical potential /i and to a 
magnetic field B : 


H^ b = H 0 -^N-2BS z , (1.3) 

with the conserved electron number N and spin S z operators defined as 



i =1 i =1 


Lieb and Wu showed that this system is integrable and that it can be solved by means of the 
Bethe Ansatz (BA) method [12]. The spectrum of the Hamiltonian (1.3) is characterized 
by two sets of complex quantum numbers (Bethe roots): the charge momenta {kj}jL 1 
and spin rapidities , {A J }^£ 1 . N is the total number of electrons and M is the number of 
down-spin electrons, with 2 M < N < L. The Bethe roots are solutions of the Lieb-Wu BA 
equations: 



e 


ikjL 



sin (kj) — ?'u 
sin (kj) + iu 


sin (kj) - iu \ _ tt ( 
sin(fcj) + iu/ ^ \\i 

m^l 


X m — 2iu 
X m + 2iu 




(1.5) 

( 1 . 6 ) 
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and the corresponding energy is given by 

N 

E = -2 ^ cos(fcj) + u (L - 2AT) - fiN - B (N - 2M). (1.7) 

1=1 

The system at finite temperature can be studied using the Thermodynamic Bethe 
Ansatz (TBA) method [15, 13]. The TBA equations for the Hubbard model, an infinite 
set of coupled nonlinear integral equations, were first derived by Takahashi starting from 
equations (1.5),(1.6) and the classification of their solutions according to the so-called 
string hypothesis [13]. An alternative approach was introduced more recently by Jiittner, 
Kliimper and Suzuki in [16], and is based on an associated integrable lattice model intro¬ 
duced by Shastry [17, 18] and the path integral approach to thermodynamics, also called 
Quantum Transfer Matrix (QTM) method [19, 20, 21, 22], One of the main results of [16] is 
a much simpler set of only three coupled nonlinear integral equations (NLIEs) of Kliimper- 
Batchelor-Pearce-Destri-DeVega type [23, 24, 25, 26]. Although the two approaches lead 
to equivalent expressions for the Gibbs free energy, a direct link between the two sets of 
nonlinear integral equations was until now not found for the Hubbard model. For the 
discussion of the equivalence between QTM and TBA in other models, see [27]. 

In this paper we shall reconsider the Hubbard model TBA inspired by some notable 
recent results obtained studying the TBA equations for A7=4 Super Yang-Mills (SYM). 
The TBA was introduced in the AdS/CFT context [28, 29, 30] to overcome the so-called 
wrapping problem [31] affecting the Beisert-Staudacher equations. Recently, this very com¬ 
plicated set of TBA equations was recast into the greatly simplified form of a nonlinear 
matrix Riemann-Hilbert problem: the Quantum Spectral Curve or P//-system [32, 33]. 
This new formulation, contrary to the TBA, treats the full spectrum on an equal foot¬ 
ing. It has already led to an impressive number of perturbative and nonperturbative 
results [34, 35, 36]. A similar reduction was also recently obtained [37] in the case of the 
A7=6 Chern-Simons theory, and made possible the determination of the so-called slope 
and interpolating functions, both nontrivial nonperturbative quantities [38]. The role of 
the Quantum Spectral Curve approach in the general integrable model framework is not 
fully understood and one of the purposes of the present work is to investigate whether a 
similar structure arises also in the context of the Hubbard model. 

One of the main achievements of this paper is the reduction of the TBA equations 
derived by Takahashi to a closed system of functional relations, involving only four functions 
P l{z), P2(z), (a = +, —), entire on a two-sheeted Riemann surface and characterised by 


a specific asymptotics. The equations are 


P^)P-(*) 

-P+(z)P^(z) = -2 sinh(20(z)), 

(1.8) 

p^QOP-W 

— Pq-(;z)P](-(z) = 2 sinh(2(/>(z)), 

(1.9) 


Tj ;1 (z) = e ~<Kz+™)+<Kz-™) ^(z), 

(1.10) 

where T<f>{z) = — l/z 2 

, and 


Tf, 1 W = P^(« + *u) 

P“ (z — iu) — P“ (z + iu) P+(z — iu), a = H, V. 

(1.11) 
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In (1.8) and (1.9), P" denotes the second-sheet evaluation of P“. The Gibbs free energy 
/ is contained in the large-z asymptotics, for instance 

inT^iCz)-^(/ + H + u). (1.12) 

We also found an alternative formulation, based on the fact that the P functions satisfy 
the functional relation 


P£(*)P&(z) = F ab (z + iu)e*« + F ab (z - iu)e~^ z \ (1.13) 


where F a j,(z) (a, b = +, —) are entire functions analytic on the whole complex plane. The 
set of equations (1.13) and (1.8),(1.9) is equivalent to (1.8)-(1.10), and may be seen as the 
analogue of the AdS/CFT P//-system. 

Furthermore, these relations imply that the zeros of the P and F functions are con¬ 
strained by a set of exact Bethe Ansatz equations. More precisely, setting Q+(z ) = 
e iBz/ u P H_( z )pH_( z ) and Q + _( z ) = e f(B-A)*/u F+ _( z ), we have 


B-i .i Q+-(«i + *u) 
Q+-(si - /u) 
e _ 2 /i Q+- {w a + 2/u) 
Q+- («j tt - 2*u) 


_ e e(si)/T, 


at Q+(si) = 0, 


Q+(w a + /u) 
Q+(w a - *u) ’ 


at Q. |— (w a ) = 0. 


(1.14) 

(1.15) 


where fi = n/T , B = B/T and e(z) = —2 T<j>(z) = ‘lizyjx — \fz 2 . Notice that the solutions 
of these equations have an infinite number of Bethe roots. In fact, (1.14), (1.15) coincide 
with the infinite Trotter number limit of the exact Bethe Ansatz diagonalising the Quantum 
Transfer Matrix. Since the latter equations are the starting point for the derivation of the 
NLIEs of [16], our analysis provides the missing link between the two different approaches 
to the Hubbard model thermodynamics. As mentioned above, this was a longstanding 
problem (see, for example, the discussion at the end of Chapter 13 of [39]). 

In this paper we also present a preliminary study of the numerical solution method for 
the ground state. At least for sufficiently high values of the temperature, we expect that 
the ground state is singled out by the requirement that all the roots Si are located on the 
second sheet. In the region of validity of this assumption, our functional relations give rise 
to a simple set of nonlinear integral equations, easy to solve numerically. We have explicitly 
checked the agreement of our results against the TBA predictions for many points in the 
region \B\ < 1, ]//] < 1, 0 < u < 2 and T > 1. However, the study of [16] indicates that 
there is a flow of roots from the second to the first sheet as T is decreased. It should be 
possible to generalise our numerical method to this parameter region, as well as to excited 
states, but we leave this problem for future studies. 

As a last remark, we point out that this infinite Bethe Ansatz is reminiscent of the 
one discovered in the AdS/CFT context in [40, 41] (however, in the latter case a relation 
with a lattice construction is not known). 

The rest of the paper is organised as follows. Section 2 contains the TBA equations 
written in a slightly modified form that highlights some of the symmetries of the model 
and can be implemented numerically without range restrictions on the parameters // and 
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B. Section 3 contains the Y-system [42] and the basic discontinuity relations [43] that 
allow its extension to arbitrary branches of the Riemann surface on which the Y functions, 
solution of the TBA equations, are defined. In Section 4, taking hints from the recent 
results in AdS/CFT and motivated by the symmetry of the TBA equations, the Y’s are 
parametrised in terms of T functions in two alternative gauges (T v and T H ). The simple 
analytic properties in these two gauges motivate a further reparametrisation of the T’s as 
2x2 determinants of more elementary objects: the P functions. Moreover, a resolvent-type 
representation for the P’s allows to express all the relevant quantities in terms of a pair of 
densities p w and p H fulfilling a new set of NLIEs defined on a finite support (described in 
Section 7). Equation (1.13) and the exact Bethe Ansatz equations are derived in Section 
5, where the connection with the results of [16] is also discussed. The free-fermion limit is 
briefly discussed in Section 6. Section 7 contains the new NLIEs together with preliminary 
numerics and a brief description of the numerical technique adopted. Section 8 describes 
the formal transformation relating the thermodynamic equations to the finite-size ones 
describing the Hubbard chain with twisted boundary conditions. The more technical parts 
of the analysis are confined to four Appendices. A proof of the special analytic properties 
of the Y functions on the Riemann section with only short cuts (the “magic” sheet of [44]) 
is reported in Appendix A. The study of the monodromy properties of the P functions 
and the proof that they live on a two-sheet Riemann surface is given in Appendix B. The 
derivation of the factors connecting the T v to the T H gauge, the proof of equation (1.10) 
and the derivation of the novel set of NLIEs are the main results of Appendix C. Finally, 
Appendix D contains a dictionary linking this work to the paper [16]. 


2. Thermodynamic Bethe Ansatz equations 


In the notation of Takahashi, the solutions of the TBA equations are 


{Vn(z),v'n(z),C(k)} , n<EN + , 


( 2 . 1 ) 


where the physical range for the arguments is z G M and k G [— 7T, 7t] . It is very convenient 
to reparametrise the variable k as z = sin(fc), and introduce the two-indexed functions 
Y mjTl as follows: 


Y 1>n (z) 
Yi, i(z) 

Y 2 , 2 (Z) 


Vn-l(z) 

VC (*0, 

C(k), 


Y n ,i(z) = 1 /v' n -i(z), n G N + , 

7T 7T" 


for k G 
for k G 


2 ’ 2 J ’ 


7T" 

u 

"7T 

— 7T,- 

—, 7T 

L ’ 2 J 


L 2 J 


( 2 . 2 ) 

(2.3) 

(2.4) 


The structure of the TBA equations is thus codified on the L-shaped diagram represented 
in Figure 1, with every node of the diagram associated to one of the unknown Y functions 
of the TBA. As a consequence of the change of variable z = sin(fc), the functions Y\^(z) 
and 1/Y 2 ^ 2 {z) have a branch cut of square root type on the real axis for z G (—1,1), and 
are in fact two branches of the same function Q(z)\ denoting with a tilde the analytic 
continuation around one of the branch points z = ±1, we have Y\ t \{z) = 1/Y 2 p(z). 
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Figure 1: The diagram on which the Y-system 
(3.4) is defined. 
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Figure 2: The extended diagram carrying the in¬ 
dices of the T functions. 


The TBA equations consistent with the Hamiltonian defined in (1.3) can be written 
as 


\nY n+l) i(z) 




1 + Y 2 , 2 {v) \ 
1 + 1 /Y hl (v)) 


In Z n+1 {fi), 


OO 

a n (v — z) dv — In 

m= 1 


/ 1 + Y m+ iq \ 

\1 + 1/^m+i (A) / 


* a 


m,n\Z) 

(2.5) 


lnYi in+ i(+) 



/ l + l/l 2 , 2 (u) \ 

V i + n,fW ) 


+ lnZ„ + i(H), 


OO 

a n {v — z) dv + In 

m =1 


/ 1 + l/Tgrn+i \ 

\1 + l/^m+l {B) J 


* a 


m,n\Z) 

( 2 . 6 ) 


OO / 

In Y 2 i 2 (z) = — In yi ; i(2:) + 2 s{z)/T = ^ ( In 

m=l \ 

+ In Z(B, jj.) + e{z)/T — 2 u/T, 


1 + ^m+1,1 A / 1 + l/l l,?n+l \\ 

1 + 1 /z m+1 {{i)) - 11 \l + l/Z m+1 (B)) ) 


* a n {z) 


(2.7) 


where we have defined 


a n (z) = 




B = B/T, /} = /r/T, e(z) = 2iz-\/l — 1/ z 2 = —2 Tcfi(z), 

n 

mi 


7r(n 2 u 2 + z 2 ) 1 


,n(^-) — ( (®m+n—2j(^ : ) + n+'2j (-2-) ) ; 

3= 1 


cosh(/i) 

cosh(H) 


, Z n (x ) — 


sinh 2 (nx) 


sinh 2 (x) 


- 1. 


( 2 . 8 ) 

(2.9) 

( 2 . 10 ) 


In (2.5)-(2.7), the convention for the convolutions is a * b(z) = f R dv a{v)b{v — z). Further¬ 
more, here and in the rest of the paper we are implicitly assuming that the integrals over 
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the interval z G (—1,1) run slightly above the branch cut. 

Notice that the TBA equations (2.5)-(2.7) are written in a slightly different form as com¬ 
pared to [13]. For the equations of [13], the standard method of iterative solution [45] 
requires sign restrictions on B and //, while here this can be avoided since the symmetries 
B -H- —B and // — n of the ground state solution are explicit. A further symmetry 

corresponds to the exchange of the two wings 


Fa. 6 G)' 1 /Ff,,a) B /X, 


( 2 . 11 ) 


up to a change of sign for all the driving terms (T <->■ —T, cj) —(f)). 

There are many equivalent expressions for the Gibbs free energy, for example: 

1 f 1 

f/T + p, = u/T - ln(2 cosh(/t)) - — I ln((l + Yi,i(z))(l + l/Y 2 ^{,z))) k (z) dz 

< 2 - 12 > 


where 


k(z) = arcsin(z) = iln(—x(z)), x(z) = iz + izyjl — l/z 2 , k^z) = —k (z) = 


dz 




k n (z) = k (z + mu) — k (z — mu), 


(2.13) 


so that k], (z) = 7 1 -|— 7 =====. The free energy can equivalently be rewritten 

yl-(z+inu ) 2 yl-(z-mu ) 2 

using only the Y functions of the horizontal part of the diagram of Figure 1, as 


f/T + fi = -u/T - ln(2 cosh(T)) - — / ln((l + 1/Yi ; i(z))(l + Y 2 , 2 (z))) k'(z) 


2 tt 


'-1 


— f; f In f 1 + V y l.n + l(f) > ] k;(z) 
Vi + i/z n+ 1 (B)y 


(2.14) 


The equivalence between (2.12) and (2.14) can be checked by using the TBA equation 
(2.7), and reflects the symmetry (2.11) and the properties [39]: 


f{-B,n,T) = f(B,n,T), f(B,-ii,T) = f(B, /q T) + 2/i. 


(2.15) 


Equations describing excited branches of the free energy - which control the correlation 
lengths among local operators at finite temperature - can in principle be obtained by 
analytic continuation [46], see also [47, 48, 49, 50]. For the Hubbard model, this has been 
accomplished only for a few states [51, 52, 53]. The TBA equations appear not to be the 
optimal tool for this analysis, and the results of [51, 52, 53] were obtained adopting the 
QTM method. It would be important to have a more complete understanding of the free 
energy spectrum, and one of the aims of the present work is to propose an alternative 
formulation which seems to have some advantages over the existing approaches. 
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-l+i2u 


l+i2u 


g 


-1-12 u f i 1 —i 2 ii 



Figure 3: The “magic” sheet Figure 4: The “mirror” sheet, on which the Y- 

system and T-system are defined. 


3. The Y-system and the discontinuity relations 


In this Section we shall describe the essential analytic properties of the system, necessary 
for the simplifications to be . An important complication comes from the fact that the Y 
functions have square root branch points at certain positions in the complex rapidity plane. 
In particular, it can be seen analysing the TBA equations that Y\ n and Y n \ are analytic 
in a strip |Im(z)| < (n — l)u, but have square root branch points at z = ±1 + i(n — l)u, 
z = ±1 — i(n — l)u, and possibly at other positions further from the real axis. 

To describe this structure, let us setup some useful notation. We adopt the convention that, 
for complex z, the Y’s are evaluated on a Riemann sheet with short cuts, the “magic” sheet 
(see Figure 3), and use the notation g(z) for the analytic continuation of a function g(z) 
around either 1 of the branch points z = ±1. Since we reserve this notation to branch points 
of square root type, we always have the property g = g. Another notation that we shall 
use often is g^ +n \z) = g(z + mu), where the shifts are evaluated on the magic sheet. 

From (2.5)-(2.7) with n = 1, one can derive the following discontinuity relations [43, 54] 
for the monodromies around some of the branch points closest to the real axis 


JMa> 


1 /Y 2t2 (z), Y 1A (z)Y 2t2 (z] 

1 + 1/Yi,i(z) 

1 + Y 2 , 2 (z) 

1 + 1 /Y 2 , 2 (*) 

1 + Yi,i(z) 


_ e -40(^) 


(3.1) 

(3.2) 

(3.3) 


Furthermore, a crucial property implied by the TBA equations is that the Y’s are solutions 


x It can be checked that the result of the analytic continuation around two branch points symmetric with 
respect to the imaginary axis is always the same. 
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of the Y-system [42] 


Y m ,n(z + *U) Y mn (z - ill) 


(i + Y m)n+ i(z)){\ + y m ,„— i(z)) 

(1 + l/^m+l,n(^))(l + ^/Y m -\ y n(z)) 


(3.4) 


where (m, ra) E {(1, fc), k E N + }U{(fc, 1), k E N + } with boundary conditions Yfc,o = l/>o,fc = 
Yk+ 2,2 = l/Y 2 ,fc +2 = 0 for k E N + (see Figure 1). Notice that there is no independent equa¬ 
tion centered at the node (2, 2). Indeed, the function Y 2,2 is simply related to Ygi through 
equations (3.1). In (3.4), we have used the notation Y to emphasise that the functional 
relations (3.4) are valid on the specific Riemann section shown in Figure 4, the “mirror” 
sheet, where all the branch cuts are traced as semi-infinite lines parallel to the real axis, 
and do not cross the strip |Re(z)| < 1. This gives a precise prescription on how to evaluate 
the shifted values appearing on the lhs of (3.4). Adapting the arguments of [43, 54, 55], it 
can be proved that the functional relations (3.1)-(3.4), together with the asymptotics 


f l + l/Y hn+1 (z) \ 

V 1 + 1 /Z n+ i{B) ) 


In 12 , 2 ( 2 ) 



/ 1 + Fn+1,1 (2) \ 

\1 + 1/2W l(A)/ 



In Z(B, p) + e(z)/T - 2u/T + O 



n E N + , (3.5) 
(3.6) 


and the assumption that the Y’s have no zeros or poles in the strip 1 1111 ( 2 :) | < u, are fully 
equivalent to the TBA. 

One might expect that each Y function should display, on a generic Riemann section, an 
infinite ladder of further square root branch points at steps of 2 iu, replicated from the 
ones closer to the real axis by the Y-system (3.4) [43]. However, this is not the case. We 
shall indeed prove that each of the Y functions has only a finite number of branch points 
on any sheet and that the number of Riemann sheets is actually finite. The fact that 
all functions appearing in this problem are defined on a finite genus Riemann surface is 
perhaps obvious from the perspective of the exact Bethe Ansatz and the Quantum Transfer 
Matrix construction of [16]. However, starting from the TBA equations this property is 
much harder to prove. To establish this result, we adopt the following strategy. First, we 
show that the Y functions have at most four branch cuts on the magic sheet: 


• Y\ >n {z) and Y n> i(z) for n > 2 have only four branch cuts at 


z E (—1,1) ± iu(n — 1), z E (—1,1) ± iu(n + 1), 


• 14 . 1 ( 2 ) and > 2 , 2 ( 2 ) have only three branch cuts at 

z E (—1,1), (-l,l)±2iu. 

This result is established in Appendix A 2 . Secondly, it is shown in Appendix B that no 
further branch cuts can appear on the other sheets. A direct numerical solution of the TBA 

2 It is important to notice that the proposed proof can be straightforwardly adapted to the known 
AdS/CFT cases, this gives an alternative and perhaps more transparent way to understand the nice prop¬ 
erties discovered in [44] for the AdSs/CFT 4 Y functions on the magic sheet. 
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Re(z) 


Figure 5: A contour plot for '$ , 2 ~i i z ) i n the complex 
z-plane confirms the presence of only four short cuts 
for Y 2 ,i(z) on the magic sheet, (u = 1, T = 0.5, 
B = 0.4, fj, = 0.1) 



-5 -4 -3 -2 -1 0 1 2 3 4 5 


Re(z) 

Figure 6: A contour plot for in the complex 

z-plane confirms the presence of only four short cuts 
for Yi, 2 (a) on the magic sheet, (u = 1, T — 0.5, 
B = 0.4, /r = 0.1) 


equations highlights, unequivocally, the presence of at most four short cuts for the functions 
Y\ n and Y n i , both on the magic and mirror sections. For the magic sheet, contour plots 
for the functions and with 



|l + (Y a , b (z)) ±1 | 
l+\l + (Y a ^z))^ 


(3.7) 


are displayed in Figures 5 and 6. 


4. The T-system 


Following the strategy which allowed the simplification of the AdS/CFT TBA, we shall 
now make a chain of simplifications and reduce relations (3.1)-(3.4) to the finite set of 
constraints (1.8)-(1.10). The first step is to parametrise the Y functions as 


_ T ntS+ i(z)T ntS -i(z) 
T n +i,s(z)T n ^x tS (z) ’ 


n,s£ N + , 


(4.1) 


where the T’s {T n s with n,s G N) are defined on the extended lattice displayed in Figure 
2 and satisfy, on the mirror section, a discrete Hirota equation (the T-system) 


T n ,s(z + iu)f, hS (z - iu) = f n+ i, s (z)f n _i tS (z) + T ntS _i(z)T nt s+i(z), (4.2) 

with boundary conditions = T_= Tk+ 2,3 = Tz,k +2 = 0 for k G N + . No¬ 

tice that the Y-system (3.4) is automatically fulfilled as a consequence of (4.1) and (4.2). 
However, the parametrisation (4.1) is not one-to-one. Different solutions of the T-system, 
corresponding to the same set of Y’s, are linked by a gauge transformation [56]. A clever 
gauge choice can greatly simplify the problem. 
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4.1 Vertical and horizontal gauges 

We shall now introduce the two gauges T H and T v , each defined by a set of simple condi¬ 
tions on one of the wings, horizontal (H) or vertical (V), of the diagram in Figure 2. We 
require that they fulfill the following properties, compatible with the cut structure of the 
Y functions, summarised in Section 3: 


and 


Tl, n W,(»>l) 

= 1 ; 

T v „ a (z),(n>l) 

T v „, 0 (^) = 1 , 


has only two branch cuts at z G (—1,1) ± inu 
Ti n (z)=r^ n] (z)r^- n] (z), (7i > 2 ) 


has only two branch cuts at z £ (— 1 , 1 ) ± inu 
T l, 2 (z) = (z) T^ [ r n] (z), (n > 2 ). 


(4.3) 

(4.4) 

(4.5) 

(4.6) 


The conditions listed above still leave a residual gauge invariance. These gauges are fixed 
uniquely by imposing the following extra conditions: 

T? i 0 (z) = e~ 2 ^ z \ T^Os) = e 2 ^ z \ (4.7) 


It is shown in Appendix B that this choice is indeed possible. Given these definitions, it is 
far from obvious that also the functions T(] m with n > m and m with m > n have such 
a simple cut structure. However, one of our main results, is that these two gauges are in 
fact connected by an elementary transformation: 


n, a ( z ) 

w 

n, a ( z ) 


(0) — 0[+ S l (2)+(/>[ (z) 

(z) ~ 


s G N + , 


1 

%Xz) 

T” 2 («) 


= g—^ En=o sgn(—s+l+2n) ^~ s+1+2n ] (z), g £ 

_ e -<fi s+1 \z)+<f>^ s - l \z)+^ s - 1 \z)-<f>^ 1 - s Xz) s __ 2 3 


(4.8) 

(4.9) 
(4.10) 


To prove (4.8)-(4.10), it is necessary to parametrise (4.3)-(4.6) with more elementary build¬ 
ing blocks, the P functions introduced in Section 4.2 and the resolvent parametrisation of 
Section 4.3. The technical details of the derivation are described in Appendix C. 

In the matching condition (4.8)-(4.10), together with the requirement that the T’s have 
no poles, is hidden the full content of the original TBA equations, and its generalisation 
to excited states. Furthermore, we will see in Section 5 that the T functions have a clear 
interpretation: they are directly related to the eigenvalues of the (fused) Quantum Transfer 
Matrix in the infinite Trotter number limit [16]. 


4.2 The P functions 


The next natural step is to parametrise the T functions defined in (4.3)-(4.7), and possessing 
only two cuts, in terms of objects having only one branch cut running along (—1,1). We 
introduce a pair of P a (a = +. —) functions for each of the wings, and write 


T l 8 (z) = P+ [+s] (z)Pl hs] (z) - P+ hs] (z)P^ [+sl (z), (4.11) 

T^ ( 1 (z) = P+ [+s 1 (z)p!_ [ " s 1 (,z) - P^~ s] (z)P v _} +s] (z), (4.12) 
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with s 6 N + . Thanks to relations of Phicker’s type among determinants, this parametri- 
sation is automatically consistent with the T-system in the corresponding H or V wing of 
the diagram. This “quantum Wronskian” construction [57, 56] suggests that the P’s are 
related to the Q functions describing the spectrum of the thermodynamic problem [16]; 
this relation is made precise in Section 5. 

To fix completely these functions, we have to specify their asymptotics and to constrain 
further their analytic properties. As suggested by the asymptotics of the Y functions (3.5), 
we demand that 3 , for large z, 

P£(z)- . Ch e ^ §z / u , P^(z) ~ Cv e^A*. (4.13) 

y 2 sinh(R) V 2 smh (A) 

The particular normalization of the prefactors C H and C v in (4.13) is chosen for future 
convenience. As shown in Appendix B, the four P functions live on a two-sheet Riemann 
surface. This is an important simplification as compared to the AdS/CFT case studied in 
[32], where the analogous P functions have additional infinitely many branch cuts, starting 
from their second sheet 1 . As we discuss in Appendix C, equation (4.7) gives the constraint 

P+(z)P!!_(z) - P+(z)P!!.(z) = -2 sinh(2(4.14) 
P+(z)P^(z) - P+{z)P v _(z) = 2 sinh(2 <j>{z)). (4.15) 

Equations (4.14)-(4.15), together with the relation giving the matching of the two wings: 

Tj ; 1 (z) = e-^ I+11 ( z ) + ^ I_11 ( 2 ) T^z), (4.16) 

the asymptotics (4.13) and information on the number of zeros on the first Riemann sheet, 
form a closed set of conditions for the P functions. In Section 7, we will show how to solve 
these equations numerically. Notice that the constants C H and C v are fixed by the solution 
to this system, and contain the free energy / through the relation: 

f/T + fi = -21nC H -u/T = —21nC v + u/T. (4.17) 

We will prove (4.17) in Section 4.4, with the aid of a very convenient parametrisation in 
terms of resolvents, inspired by [44, 41]. 

4.3 Resolvent parametrisation 

In the following, we shall assume that the P functions have no zeros on the first sheet. 
We expect that this singles out the ground state solution in a large parameter region for 
sufficiently high temperatures. The numerical solution of the functional relations presented 
in Section 7 reveals the existence of an infinite number of zeros on the second sheet. 

3 In analogy with the AdS/CFT case [33], for excited states we expect that this asymptotic behaviour 
may be modified by power-like prefactors encoding the quantum numbers. 

4 In AdS/CFT, the structure associated to the presence of a ladder of infinitely many branch cuts is 
encoded in a periodic matrix fi l: j. In the current setup, the introduction of a periodic matrix can be 
avoided. 
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Starting from the ground state, we expect that excited states can be obtained by analytic 
continuation in B and / 1 , leading to the migration of a finite number of zeros to the first 
sheet and to simple modifications in the equations presented in this Section. We point out 
that, in a neighbourhood of T = 0, even the ground state solution should present a finite 
number of zeros on the first sheet [16], and therefore the following analysis needs to be 
modified. 

Let us introduce two resolvent functions 

«*« - h /_' 7^v iv ’ Q = H ' V ' <4 ' 18) 

where the densities p H and p v are for the moment undetermined. This is the most generic 
parametrisation of a function with a single cut, no zeros and no poles on the first sheet, 
and large-z asymptotics ~ 1/z on the first sheet. Under analytic continuation through 
the cut, the resolvents transform as 

Ga(z ) = Sa(z) + Pa{z). (4.19) 

Thus, under the assumptions discussed above, we can parametrise the P’s as 

\J~2 sinh(B)P!j.(z) = h H (*) e ±(BG H (*HB* /u ), (4.20) 

\/ r 2 sinh(/t) P^(z) = h v (z) e ±(Ae v U)-f A^/ u ) ; 

where the functions h Q (a = H, V) are required to have a single cut, no zeros or poles and 
a constant leading asymptotics on the first sheet. Taking ratios of these functions we find 

, )R P"P H „ P[),P V 
e 2Bp " = ~ + - e 2 MP v = ~ + (4.21) 

P H P(j_ P v P^_ 

As a consequence of the simple monodromy properties of the P’s, also p H and p v live on a 
two-sheet Riemann surface, with 


Pa(z) = -p a (z). (4.22) 

More precisely, we require that the behaviour at the branch points is such that the functions 
Pa(z)/(z^l - l/z 2 ) are analytic in a neighbourhood of the cut z € (—1,1). Due to (4.21), 
however, the densities can have logarithmic singularities elsewhere in the complex plane, 
in correspondence to zeros of the P’s. 

Substituting (4.20) into (4.14) and (4.15), we find 

h H (z)h H (z) sinh(Rp H (2 : )) = 2 sinh(R) sinh(2(£(z)), (4.23) 

h v (2)h v (z) sinh( / u/> v (2 : )) = —2 sinh(/i) sinh(2(()(z)). (4.24) 
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PH 


~PV 




Figure 7: p H for B = 1, p, = 0.5 and T = 1, for Figure 8: — p v for B = 1, p = 0.5 and T = 1, for 
u = 0, u = 0.25, u = 0.5, u = 1 and u = 2. p H is u = 0, u = 0.25, u = 0.5, u = 1 and u = 2. —p v 
monotonically increasing with u. is monotonically decreasing with u. 


The solution of the latter equations, compatible with the asymptotics (4.13) and the ab¬ 
sence of zeros on the first sheet, is unique and given by 5 : 


In h H (z) = 

lnhv(z) = 


ZyJ\ - l/z 2 I" 1 ln ( 2sinh(2</>(t>))sinh(ff) \ dv ^ 

2vr ,/_i y sinh(l?/9 H (u)) J \/l — v 2 (v — z) ’ 

zV 1 - V~ 2 f 1 2 sinh(2^(u))sinh(/i) y du ^ 

2 vr 7_i n V sinh(/ip v (u)) / y'l - u 2 (u - z) ’ 


Thus, with the dehnitions (4.25) and (4.26), the parametrisation (4.20) automatically fulfills 
the constraints (4.14),(4.15). 

Considering the large-z asymptotics of (4.25),(4.26) we find that, on the first sheet, 


lim lnh H (z) 

z—t oo 


lim lnh v (z) 

z—> OO 


1 j" 1 1 ^2sinh(2^(u)) sinh(S)^ dv 

2 vr i y sinh(5p H (^)) / \/l — u 2 ’ 

1 f 1 1 f 2sinh(20(u)) sinh(/i)^ dv 
2 vr ,/-i n V sinh(/lpv(u)) / y/l - u 2 ’ 


(4.27) 

(4.28) 


In Section 4.4, we will connect the integrals appearing in (4.27) to the free energy and 
establish relation (4.17). 

From (4.11)-(4.12) and (4.20), one can rewrite the T functions in terms of the densities as 


TpH id+dl, C s ]-- 7-H rjpV Vd^Jvd - s JtV 

M,S — n H n H /1,S) A s,l - n v n v / s ,l 


— v,[+dv.[-s]' 


G N, 


(4.29) 


where we have borrowed the notation of [44, 41] and denoted 


T H — 

n,s ~ 


sinh \B(s + G l „ l - G l ~ 
sinh (13) 


T v — 

's,l - 


sinh 


s + Q? - Gl 


sinh(/t) 


(4.30) 


5 Notice that in (4.25), needs to be evaluated just above the branch cut, so it agrees with ^\/l — 
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Figure 9: Comparison of the numerical solution of the TBA equations (red, dashed line) with the 
new set of nonlinear integral equations (gray, continuous line), atT = l,u = l,£? = 0, /r = 0. The 
two curves are interpolations of the numerical solution. The numerical solution of the TBA was 
obtained by discretising the integrals with a step Az = 1/60 ~ 1.7 x 10” 2 , truncating the infinite 
sums after 60 terms and cutting off the integrals over the real line at z max = — 2mm = 50. The 
new NLIEs were solved as explained in Section 7, using lV trunc = 50 in (7.8). The difference of the 
curves is of order 1CF 6 , and is not visible on the scale of the plot. 


Equation (4.29) defines a gauge transformation between the T" and the T a functions, with 


Tq,s ~ 1) 7^o — 1) s G N, 

_ n~ H [+ s l •t h [~ S 1 t-v _'t'V [+ s ] -t- v [ —s ] ,2 m+ 

'2 ,s ~ '1,1 '1,1 ) 's,2 ~ '1,1 '1,1 i S fc 1^1 . 


(4.31) 

(4.32) 


In particular, the h factors in (4.29) cancel in the gauge invariant combinations (4.1). The 
gauges T" are characterised by the simple large z asymptotics 


Ti s ( z ) ~ sinh(sB)/sinh(5), T^i(z) ~ sinh(s/i)/sinh(/i), s G N + , (4.33) 


which fixes them uniquely. Finally, in these gauges relations (4.7) become 


K 0 (z) = 


sinh (Bp H (z)) 

2 sinh( 2 (/>(z)) sinh(B) 


=-2 


ToM = 


sinh(/ip v (z)) 


2 sinh( 2 (/)(^)) sinh(/t) 


d 2 ^(«) 


. (4.34) 


As discussed in Section 7, the densities can be computed by solving numerically a set of 
NLIEs, which can be derived from the constraint (4.16). Some examples of solutions are 
displayed in Figures 7 and 8 . Finally, let us remark again that the information contained in 
the pair {p H , p v } is fully equivalent to the knowledge of the solution of the TBA equations. 
In Figure 9, the TBA solution is compared with the same quantity reconstructed from the 
parametrisation (4.1),(4.30). 


4.4 The free energy 

The purpose of this Section is to derive a simple formula for the free energy in terms of 
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the densities. The starting point is equation (2.12), written in the equivalent form: 


f/T + fi = 

if/ 

2vr ^ L 


— u/T — In (2 cosh (5))- 

27T 


'-1 


In ((1 + 1/Yi,i(z))(l + Y 2}2 (z))) k '(z) dz 


n= 1 ' 


, ( r i,n +1 ( Z + iu )V!n+l ( Z ~ * U ) \ , / , x , 

ln l —— Ki > • 


(4.35) 


where we have set Tf s = T/f s sinh(.B)/(sinh(si?)), so that In T" s (z) ~ 0 at large z. Adopting 
the “telescoping” technique used in [55, 44, 40], it is now possible to remove the infinite 
sum appearing on the rhs of (4.35). Since 7) H n is analytic for |Im(z)| < (n + l)u, we can 
shift the integration contours and prove that, for n 6 N + , 

[ infej+l ] (z)T^l ] (z)) k' n {z) dz = f ln(T" n+1 {z)) (k^ +1 (z) + k! n ^{z)) dz. 

J R v 7 JR 

(4.36) 


This shows that many cancellations take place in (4.35), leading to 

1 


f/T + A = -u/T - 


2 vr 


J 1 In ((1 + l/Y hl (z))(l + Y 2 , 2 {z))) k \z) dz (4.37) 


- [ lnT" 2 (z)k , (z)dz + / hi i(z) k/ (z) dz. 

ft J -1 47r 


Using (4.1) and the fact that > 2,2 = 1/Tpi, (4.37) becomes 


1 


f/T + fi = -u/T - — J ^ in 


1 


1 I T?J +1] (z) T"t 1] (z)T"J +1] (z)T"i~ 1] (z) 


V 


Ko(z)Tf 0 (z) 


+ — / lnTf ] 1 (z)k[(z)dz. 

./to 



The convolution appearing in the first line of (4.38) can be viewed as a contour integral 
around the cut of k!(z), and, deforming this contour to a pair of infinite lines, one can 
write 

fh In (7^ 1 I+I| W7^ 1 , - 1| W7^ 1 ,+I| W7^ 1 [ - 1| W) k \z)dz 

= ([ ) ln (7?,i +1] { z )7i,l~ 1] { z j) k\z)dz. (4.39) 

\JR+iO+ JR-iO+/ v 7 

Shifting the line contours appearing in the latter expression (and pushing some of them to 
infinity), one can show that the integrals involving Tfi in equation (4.38) actually cancel 
each other out. Substituting (4.34) for 7/ H 0 , we finally arrive at the desired result: 


T f 1 /4 sinh 2 (2<)>(z)) sinh 2 (.B)\ dz 

f+^ = -n-—J ^ny J -j—. 


(4.40) 
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Starting from the alternative expression (2.14) for the free energy, working in the gauge 
T v and following the same steps, we also obtained 

T /‘inf t™h 2 (20(z)) 3 mh 2 (A) \ dz (441) 

2 tt V sinh 2 (/f/9 v (2:)) / \/l - z 2 

The comparison between (4.40), (4.41) and (4.27) finally gives the formula (4.17), showing 
that the free energy / appears in the asymptotics of the P’s. 


4.5 Energy-carrying Bethe roots 

Let us look more closely at the kernels appearing in (4.40) and (4.41). From (4.21), we 
have 

smh 2 (2 <j,{z)) = p^) P H^)pH^) P H (.) (4.42) 

sinh 2 (Bp H (*)) + + 

where we have used (4.14). The same relation with B —> p, H —> V is valid for the 
density p v . Thus, the kernels entering the energy formula have logarithmic singularities 
in correspondence to the zeros of the P functions. Should any of these zeros cross the 
integration contours in (4.40),(4.41), the energy would pick extra residue terms of the form 
iTk(zj), as is easy to verify by first performing an integration by parts 6 . This suggests 
that zeros of the P functions can be regarded as energy-carrying Bethe roots for the finite 
temperature spectrum. According to our assumptions, for the ground state with real B 
and p all the zeros lie on the second sheet, but they may cross the interval and move to the 
first sheet under analytic continuation in B and p. We will return on this point in Section 
6 below. 


5. Exact Bethe Ansatz and relation with the Quantum Transfer Matrix 

In this Section we will recast the closed set of conditions (4.14)-(4.16) into another, inter¬ 
esting form, which reveals the presence of an exact Bethe Ansatz. This second formulation 
is perhaps closer in spirit to the P/j-system of AdS/CFT [32], which was one of the main 
sources of inspiration for the present work. 


5.1 Formulation as a coupled Riemann-Hilbert problem 

We start by considering the analytic continuation of (4.16) through the cuts at zG (—1,1)± 
m, which yields 


P" + 1±2| (z) fe(i) - P" l±21 (2) P+(f) = 
p)y [±2] (z) PL (z) - p^ [±2] (z) P^ (z) 


(5.1) 


Two of the P functions in (5.1) can be eliminated using the constraints (4.14), (4.15). 
Solving (5.1) for the remaining two P’s, we find 

+ a, !>=+,-. (5.2) 

6 One could also express the ground state free energy as an infinite sum iT L), k (zj) (which should be 
regularised). 
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where F a & has a simple explicit expression: 


F ab(z) = 


z* [ 1] ( z )Pa [+11 (2)Pj [ l! (z) + e^ 1+1] Wp" [ 1] (z)P^ [+11 ( 

T?,i( 2 ) 


a,b = +,—.(5.3) 


Relation (5.3) shows that F a b(z) could - in principle - have a pair of branch cuts on the lines 
Im(,2) = ±iu. However, by shifting equation (5.2) of ±iu, and evaluating the discontinuity 
across the branch cut, we find 




disc 


r [±i] 

ab 


= disc 


pH [±1] p v [=1=1] 


= o, 


(5.4) 


which implies that F a b is analytic in the whole complex plane. The absence of branch cuts 
can also be directly verified by using the definition (5.3): the quantities F^ — f[^ vanish 
since they are proportional to the difference of the lhs and rhs of (4.16). Notice that the 
lhs of (5.2) has no poles. As a consequence, it is possible to prove that the F functions, 
appearing on the rhs, are entire on the whole complex plane. This can be also deduced 
from equation (5.9) below. 

These newly introduced quantities fulfill, together with the P’s, a closed set of func¬ 
tional relations. The fundamental set of equations is 


P+(z)P^0?) - P+(z)P!!.(z) = -2 sinh(20(z)), (5.5) 

P+(+)P^(z) - P+(z)P^ (z) = 2 sinh(2 c(5.6) 
P^)P^) = F't+e'W*! + «,(> = +,-, (5.7) 


with the requirement that the F functions have no cuts and no poles, and the P’s have 
only one cut on each of their two sheets, and no poles. These equations make explicit the 
nonlinear Riemann-Hilbert nature of the problem, since they show how the two branches 
of the P’s are connected through the entire 2x2 matrix F a ^. In this respect, they are 
reminiscent of the P//-system. It is important to underline that (5.5)-(5.7) are completely 
equivalent to (4.14)-(4.16). In particular, it is possible to show that the definition (5.3) 
and the relation (4.16) are both hidden in these equations. 

Equations (5.5)-(5.7) imply the existence of many other functional relations among the 
F’s and the P’s. For instance: 


F ++ (z)F__( 2 )-F+_(z)F_+(z) = -1 (5.8) 

F !+^(z) F [ Z^(z) + F^i ] (z) F^(z) = f£I ] (z) F[j(z) + F ^(z) F l ~%) + 2 cosh(20(z)). 

Finally, as already remarked in Section 1, the Hubbard Hamiltonian (1.1) has, for B = 
H = 0, a hidden 50(4) = 51/(2) x 517( 2)/7,2 symmetry [14]. While this symmetry was not 
very evident in the original Y-system and discontinuity relations, it appears to be nicely 
encoded in the structure of equations (5.5)-(5.7). At generic values of the magnetic field 
or chemical potential, the symmetry is broken by the boundary conditions (4.13). 
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5.2 The exact Bethe Ansatz 


The relations obtained in the previous Section contain an exact Bethe Ansatz constraining 
the position of the zeros of the P and F functions. Starting from (5.2), we immediately 
get 



(*) = ± 


1 P2(z)Pb(z) - e^«P"(z)P£(z) 

2 sinh(2(/)(z)) 


a, 6 — . 


(5.9) 


Combining these expressions to form bilinear combinations of the F a f,’s reveals a further 
set of interesting relations of quantum Wronskian type. For instance, using (5.9) and 
(4.14),(4.15), we obtain 

F^WF^'W - FH( Z )FM( 2 ) = P|j_(z)P!j_(z). (5.10) 


The rhs of (5.10) defines an entire function on the whole complex plane. It is convenient 
to isolate the exponential prefactors and set 

Q H + (z) = e^ z / u P H + (z)P H + (z), Q + -(z) = e^-^ u F+_{z). (5.11) 


The a = +, b = — case of equation (5.2) can then be rewritten as 





(5.12) 


We assume that for the ground state all zeros of the P functions are on the second sheet; 
therefore the zeros of the rhs of (5.12) coincide with the zeros of Q(j_, giving the first of 
the Bethe Ansatz equations quoted in Section 1. The second BA equation can be derived 
by taking ratios of the expressions obtained by evaluating (5.12) at zeros of Q+l. The 
resulting set of quantisation conditions is: 


e B-a Q+-{sj + m) 
Q+-(si - *u) 
e _2/i Q+- {w a + 2m) 
Q+- (w a - 2iu) 


_g-2 #(*») 


at Q+(sj) = 0, 


Q+(w a + iu) 
Q+(w a - iu) ’ 


at Q- 1 — (w a ) = 0. 


(5.13) 

(5.14) 


Very interestingly, these equations appear to be the infinite Trotter number limit of the 
BA diagonalising the Quantum Transfer Matrix of Jiittner, Kliimper and Suzuki [16]. We 
remind the reader that the QTM is defined as a discrete object acting on a IV-site quantum 
space (N is known as the Trotter number), and that the thermodynamics of the Hubbard 
model is described by the largest eigenvalue of the QTM in the limit N —»• oo. For the 
ground state at finite even values of N, the BA of [16] is 


e B-ji qijsi + iu) 
q 2 {si - *u) 
e — 2 A g 2 (w a + 2iu) 
q 2 (w a - 2iu) 


-bfr(si), i = l,...,N, 

qi(w a +iu) - 

-7-rr, a = 1,..., A/2. 

?iK - *u) 


(5.15) 

(5.16) 
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where' lim i Y_ 5l0O bft(z) = e 2< ^ and 

N N/2 

91 ( z ) = ~ s *) > Q2(z)= Yi & - w <*) ■ 

i=l a=l 


(5.17) 


The functions appearing in the BA equations (5.13),(5.14) are the continuum version of 
these polynomials. In fact, we believe that they can be factorised over their zeros as infinite 
products of the form 8 


Q+-( z ) 

Q+-(o) 



#>) 




Q H +(*) 

2 +( 0 ) 




(5.19) 


where Re(s_j) = — Re(sj), Re(ic_ a ) = —Re(rc a ). In (5.19), the zeros have been paired 
up in order to make the product convergent. This rearrangement is necessary, since the 
zeros accumulate at infinity with an evenly spaced asymptotic distribution (see Sections 6 
and 7). The regularisation (5.19), without extra Hadamard factors, is consistent with the 
infinite Trotter number limit of (5.15)-(5.17). 

In addition to (5.13), (5.14), there are naturally other equivalent sets of BA equations. 
For completeness, let us discuss the general structure. From (5.9) and (4.14), (4.15), one 
can obtain the following quantum Wronskian-type relations 


fLIWfLIW - Ffcl(*) f£](*) = ±Pi(z)Pi(z), (5.20) 

FLIWFgW - FaWFjJ(z) = ±Pl(z)Pl(z), (5.21) 

F^(z)F[rI ] (z) - F^(z)F^I ] (z) = P^(z)P^( z ) + P H + (z)P H _{z), (5.22) 

Ffi ] (A> F -+(*) - F^I ] (z)F^(z) = P\(z)P^_(z) - P^(z)P+(z), (5.23) 

which allow one to derive alternative pairs of BA equations. Generalising (5.11), we define 

Q"(z) = e i6 ^ n P H a (z)P H a (z), (5.24) 

Ql{ Z ) = <? ila l / ' 1 K(z)K(z), (5.25) 

Q ab (z) = e UBa+Mz/u Fa6 ( z ) > a, b = +, -, (5.26) 


where B± = ±B, (i± = ±£i. We expect that all the Q ab and Q a functions thus introduced 
admit a factorisation of the form (5.19), each with a different set of zeros. The systems of 

7 See Appendix D for more details. 

8 We expect that the P functions can be factorised in the following form (cf [41]), 

OO 

(2sinh(B))l e ± ^ z/u P£( 2 ) = C H (* + x(z))) (l + 1 /(z^ x(z))) , (5.18) 

i= 1 

(H +) (H +) 

where z\ ' ; = x( s\ '), and the Zhukovsky map x(z) is dehned in (2.13). Again, the pairing of zeros is 
to guarantee the convergence of the product. 
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BA equations fulfilled by these functions can be easily obtained from (5.13) by symmetry. 
We can schematically summarise these symmetries as 


Q-(z', B,jj) = Q v _(z;B,fi) = Q v + (z;B,-p), (5.27) 

Q-a(z;B,p) = Q +a (z;-B,fi), Q a -(z; B, fi) = Q a+ (z; B,-ft), (5.28) 

with a = +, —. In addition, we have the substitution rule 

Qa ^ Q'ai fi ^ fi = — fii (! = +,—• (5.29) 

Finally, let us point out that the definitions (5.24), (5.25) can be inverted as 

lnP a(z) = -\B a z/u+ [ \nQ"(v)-—J^- --, (5.30) 

lnP a(z) = ~l^az/u + v/l ~." [ In gV(u)—==J---, (5.31) 

2 2m J_ x V1 -v 2 (v- z) 

(a = +, —), leading to the following equivalent formulae for the free energy: 

/ + n = -u - - [ ln(2sinh(S) Q±(u)) dh = u - - f ln(2sinh(/}) Q±(v)) dv 

TT J -1 V 1 ~ V 2 7T J_i Vl ~ V 2 

(5.32) 


5.3 Relation with the Quantum Transfer Matrix 

By comparison with the results of [16] we can uncover the physical meaning of T" 1; T]) 1 and 
show that these functions are simply related to the eigenvalues of the Quantum Transfer 
Matrix. Starting from (5.2), we get 


^ e -2 V +%) + 
( e ^ l ~ 1] {z) + 


F ab(z) ) 

F ab (z) ) 


td h . 

0 -^ +1 \z) * a 


F ab(z) 


Fafe(^) 


v[-l]. 


a, 6 — . 


(5.33) 

(5.34) 


Multiplying (5.33) by ^/P ( H J +1 ^ and (5.34) by P^ +1 '/P^ ^ and adding them, we 
find, after using (5.3) 


>H [ 1] 


(*) 


p« [+1, w 


e -2</- [+1] (z) _|_ 


p[+ 2 l 

ab 


(*) 


>H [+1] 


F ab(z) 


+ 


(*) 


Pa' 1 ] (Z) 


= 2</>[ _1I ( z ) + F^ '(z) 


F ab(z) 


Pn [ 1] (z)P v b l+1] (z)e-^ l+1] ^+ Pa [+1] (z)Pl [ I] (z)e </>1 11 W 

F ab(z) 


= T?,iW. 


= (5.35) 


(5.36) 


The combination in (5.35) agrees with the form of the Quantum Transfer Matrix eigenvalues 
(see equation (15) in [16]), namely 


(A)ref [16] ^ TTgl- 


(5.37) 
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The precise details of this identification are given in Appendix D. Notice that, as a con¬ 
sequence of the absence of poles for the P functions, this quantity does not have poles on 
any sheet, and that this pole-free condition gives precisely the Bethe Ansatz. Furthermore, 
setting = Q", we can write 


TpH P>H [+1] pH 

11 1.1 r a r a 


[- 1 ] 


Qa [+1] 


e -20 [ + 1] + 


^[+ 2 1 ' 
ab 


ab 


+ 


Qa [_1] 


e 2^~ 1] + 


P [-2]' 

ab 


ab 


.(5.38) 


In the quantity in the square brackets we recognise the infinite Trotter number limit of 
the “auxiliary” transfer matrix eigenvalues A aux introduced in [16], equation (25). Setting 
a = +, the identification is 


(A aux ) ref j 16 j <r. 


TVi 


p H [+i] p H [ i] 


(5.39) 


6. The free fermion limit 


The exact solution of the TBA equations at u = 0 was found already by Takahashi in 
[13]. It is interesting to recover this result starting from relations (4.14)-(4.16). In the 
limit u —>• 0 + , the shifts in the P functions shrink to zero and, for Re(z) G ( —1,1), the P 


functions collapse to their values above/below the cut. Therefore, 

Ti ;1 ~ e s P+ PH. — e~ B ~P^_ P[]., (6.1) 

~ e p ¥\ PL - e" A PL P^, (6.2) 

and equation (4.16) becomes 

e BpH_ pH_ _ e -BpH_ pH_ = ^ e Apv_ pv_ _ e -M P v_ pv_^) e - 2 4 >' (6.3) 

A second constraint is simply obtained by continuing (6.3) to the second branch: 

e hpH_ p h_ _ e -BpH_ p h_ = ^ ; /ipv_ p\^ _ e -Apv_ pv_^ e 2 L (6.4) 

Solving (6.3),(6.4) with the aid of (4.14),(4.15), we find 

sinh(R) P+(z)PL(z) = 2cosh(0(z) — (B + fi)/ 2) cosh(^>(z) — (B — jj)/ 2), (6.5) 

sinh(/i) P+(2 )PL(;z) = 2cosh(^>(^) + (B + fi)/ 2) cosh(</>(z) — (B — fi)/ 2). (6.6) 


We can now compute the densities using (4.21). For the horizontal wing the result is 

e 2 b Ph {z) = P+(~) P -Q) = cosh(<ft(z) + (B + AVg) cosh(0(g) + {B - fi)/ 2) 
PL(z)P+(z) cosh(<^(z) — (B + /t)/2) cosh(^(z) — (B — p,)/2) ’ 

and, denoting the rhs of this equality as e 2BpH ^ = R(B,£i;z), the density characterising 
the vertical wing is given by 

e -2£pv« = R^B\z). (6.8) 
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Figure 10: Position of the zeros of P 1 / in the Figure 11: Position of the zeros of P^_ in the 


+ 


complex z plane, for the free fermion solution at complex z plane, for the free fermion solution at 


B = 1 and p = 0.5 and T = 1. 


B = 1 and p = 0.5 and T = 1. 


The kernel appearing in the free energy formula (4.40) reduces to 

f sinh 2 (2())(z)) sinh 2 (i?)\ ^ / sinh 2 ( 24 >(z)) sinh 2 (/i)\ 

\ smh 2 (Bp H (z)) J V sinh 2 (pp y (z)) ) 

= — ^2 In ^cosh(</>(z) + aiB/2 + a 2 p/2)^j , (6.9) 

<71,02 

with <J\,a 2 = ±1, giving the well-known result for the Gibbs free energy at u = 0. 

The pattern of Bethe roots displayed by the free fermion solution is interesting, as the 
numerical solution for u > 0 (see Section 7) suggests that the zeros are smoothly deformed 
away from their positions at u = 0. Each of the P functions has two infinite strings of 
zeros on the second sheet, corresponding to the two factors on the rhs of (6.5). Denoting 
as z n (B , p) the solution of the equation 

2 4>{z n (B, p)) + B + p = i(2n + l)7r, (6.10) 

and defining 

A(B, p) = | z n {B, p), n E zj , (6.11) 

the distribution of zeros is summarised in Table 1, and illustrated in Figures 10 and 11 
for two of the P’s at if = 1 and p = 1/2. Zeros belonging to A(B,p) accumulate at 
infinity along the line Im(z) = —B/2 — p/2, and their asymptotic spacing is ttT. As can 
be obtained from equations (5.20)-(5.23) in the u —> 0 limit, in the free fermion case the 
zeros of the F functions are a subset of the zeros of the P’s (see Table 1). 

Let us make a short comment on the analytic continuation mechanism governing the 
transition to excited states. As already anticipated in Section 4.5, the free energy acquires 
an extra residue whenever one of the energy-carrying Bethe roots crosses the integration 
contour. The movement of these zeros can be driven by analytic continuation in B and 
p. In general, such a crossing does not correspond to a branch point in the domain of 
the parameters B or p, since the contour can be deformed to avoid the contact with the 
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Function 

Zeros 

K(z) 

£ G A(B a , p) U A(B a ,-fi) 

KM 

Z G A(B, Pa) u A(-B, Pa) 

F ab(z) 

z G A(B a , -p b ) 


Table 1: Distribution of the Bethe roots among different P and F functions, with a,b = , 

and B± = ±B , p± = ±p. Zeros of P’s live on the second sheet. 


wandering zero. Genuine branch points correspond to the so-called pinching phenomenon, 
when a pair of zeros collide on the contour from opposite sides [46]. In the free fermion 
case, this happens at values of B and p given by: 

X B /X^ e ±2 / T = -1, X B A m e ±2 / T = -1, (6.12) 

where X B = e B / T , A M = are the fugacities, corresponding to a pair of zeros of 

sinh 2 (i?)/sinh 2 (ilpH) = sinh 2 (/i)/sinh 2 (/i/)v) pinching the contour of integration at the 
origin in the z-plane. Analytic continuation around one of these points causes the transi¬ 
tion to an excited branch of the free energy, and of the P functions. Finally, it is interesting 
to notice that the branch points (6.12) mark the boundaries of the four phases of the sys¬ 
tem at T = 0. In the interacting regime u > 0 at T = 0, the phase diagram includes a 
fifth phase describing the Mott insulator behaviour [12] (see also Chapter 6 of [39]). It 
would be very interesting to investigate the branching structure of the free energy at finite 
temperatures and coupling, and link it to the phase diagram of the Hubbard model. We 
expect this to be possible with the numerical method described in Section 7 and plan to 
come back to these questions in the future. 


7. Numerical solution 


In Section 4.3, we have parametrised the P functions appearing in the problem in terms 
of the densities p H and p w . They can be computed by solving a system of coupled non¬ 
linear integral equations, which determine simultaneously the densities and the function 
j entering the TBA equations. This formulation is very similar to the set of NLIEs 
proposed in [44] for the AdSs/CFT 4 spectral problem. We discuss here only the ground 
state equations. 


7.1 Nonlinear integral equations 

First, by expressing Y\\ in the T» and in the % gauge, we find 


r(z ) 


tZYA Tp-\z) 

Tj“, [+1 (2) 


= g^O) 


Tpp'(z) tY'(z) 


where 


r(z) 


( l + yi,i(z) e 4< ^\ 

1 + Y ltl {z) J’ 


(7.1) 


(7.2) 
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and the functions Tj* j , i = H,V depend on the densities through the parametrisation 
(4.30). In an iterative scheme, equations (7.1) can be used to update the values of the 
two densities starting from the knowledge of r(z). The numerical method is described in 
[44] and is reviewed below in Section 7.2. To close the system, there is a further equation 
determining Y\\ as a function of p H and p w : 


In 14 , 1 ( 2 :) 


/ 712(2) \ 

2u/T-2^) + ln^^ M j 
, ( T“! +,I (z) 

“Urw^™' 



dv 

2Tri{v — z) 


disc 


/' 7ti l+11 w'^ 1 1 ~ 11 w y 

(7.3) 


This relation is derived in Appendix C, both from the TBA equations and also purely 
from the Riemann-Hilbert formulation (4.14)-(4.16). An equivalent form of (7.3), which is 
convenient for the numerical implementation, is 


In Ti,i( 2 :) 


2u/T — 2 4>(z) + In 




dv 

2tt i{v — z) 


disc 


/ t!l +l ('Ki H (») V 



V’U) 

vT 1 ' 


( 2 ) ^'-"(2) J 


(7.4) 


where f denotes Cauchy’s principal value integral. 


7.2 The numerical method 

The system (7.1),(7.4) can be solved iteratively for the values of p H (z), p v (z) and 14 , 1 ( 2 ) 
on the interval 2 € (—1,1). One iteration step, updating the values of the densities p^ —> 
p\ k+1 \ can be represented as 

yff -^' 711 pf +1) ', A k+iy . (7.5) 

After each step, we update the solution as p) fc+1 ' — + (1 — #)/4 +11 . Tile introduction 

of the weights 6 , 1 — 9 is a common recipe used to ensure convergence. In all cases we 
considered, the scheme was stable taking 6 = \. 

In the above described procedure, the most difficult step is the solution of (7.1) for the 
densities. Let us review the method discussed in [44], concentrating on the equation for 
the horizontal wing. Using (4.30) and the basic property Q H = Q H + p H , the first equality 
in (7.1) can be written as 

, -j = sinl i(B{g [ H 2 \z) - 0 H (z) + Ph(z)/2 + 1)) sinh(.g(0 H (z) - G [ » 2 \z) + p H (z)/2 + 1)) 

sinh(R(C/i + 2 l ( 2 ) - 0 H (z) - p H (z)/2 + 1 )) sinh(R(^ H ( 2 ) - G [ h~ 2] (z) - p»{z)/2 + 1 )) 

(7.6) 
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(7.7) 


where <$ H denotes the Cauchy principal value integral 


(z) = — f 
2TTiJ-l 


*Yl dv 

z — V 


By extracting 9 p H from the rhs of (7.6), we obtain the density in terms of $ h and r. 

This equality is used to update the value of p H in the last passage of (7.5). 

The numerical evaluation of the singular integrals appearing in the NLIEs (7.1),(7.4) 
can be performed very efficiently using a Chebyshev expansion. To optimize the numerical 
method we found it convenient to discretise the densities by using a Chebyshev expansion 
of the second kind 


_ N trunc 

(C m b) = £ 

n=0 


c^U 2n (z), z€(- 1,1), a = H,V, 


(7.8) 


(where we have taken the correct parity into account) and evaluate principal value integrals 
using the properties 



dv 

7 t(v — z) 


\A - V 2 U n (v) 


~T n +i(z), 


dv 


-i ir(v — z ) y/l — v 2 


T n (v) = U n -i(z ), 


where T n and U n denote the Chebyshev polynomials of the hrst and second kind, respec¬ 
tively. To produce the data presented in this paper, we took -/V trunc = 50. In the vast 
majority of the cases we considered, less than 30 iterations were sufficient to achieve con¬ 
vergence of the coefficients centering (7.8) on the fifth digit. We observe that the error 
is approximately halved at every iteration. 

Let us also make a comment on the region of convergence. We observe that, for fixed 
values of B, p and u, the iterative scheme is convergent for sufficiently high temperatures 
- in particular a preliminary study suggests that, for arbitrary B,/i£K, the convergence 
region probably includes 0 < u < 2 and T > 1 - but breaks down below a certain threshold 
temperature. Lowering u, the breakdown temperature decreases, and this hints that, for a 
given value of T, the method should be applicable without modifications in a nonvanishing 
neighbourhood of the free fermion point. We suspect that the breakdown of the method 
for low temperature or strong coupling is related to the appearance of zeros on the first 
sheet for the ground state solution [16]. We plan to come back to this issue in the near 
future. 


7.3 Exploring the complex plane 

The numerical method we have described computes p H (z), pv(z) and Y\p(z) for — 1 < z < 1. 
Once we have a solution on the interval, we can reconstruct the behaviour of these functions 
in the complex plane. We do this in two steps. First, we use equation (7.3) to compute 
Y\ \ (z) for z on the whole first Riemann section - with cuts at z G (—1,1), z £ (—1, l)±2iu 
- from the values of the densities on the interval. From the same information, we can also 
compute Gh(z) and Gy(z) for an arbitrary complex value of z. Then, one can obtain the 

9 Notice that (7.6) can be written as a quadratic equation for tanh(£>pH/2). 
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1.0 



- 1.0 


Figure 12: Plot of the function —Re 


(ln(sinh 


in the complex z-plane, for B 


1, p = 0.5, 


T = 0.3 and u = 0.1. The positive peaks correspond to energy-carrying Bethe roots, while the 
negative peaks on the real axis are zeros of sinh(2(/>(z)) (cf. equation (4.42)). 


values of Ph{z), p^(z ) for any complex values of z by inverting (7.1). The complex zeros of 
sinh( / o H (-2 ; )) are visible in Figure 12 for the numerical solution corresponding to p = 1/2, 
B = 1, u = 0.1, showing a clear qualitative similarity with the free fermion case. 


8. Mirror equations 


In this Section we shall describe the finite-size versions of the functional relations (5.5)- 
(5.7), which form a set of equations fully equivalent to the Lieb-Wu quantisation conditions 
(1.5),(1.6) with more general twisted boundary conditions. Comparing the thermal BA 
(1.14),(1.15) and the Lieb-Wu equations (1.5),(1.6), and considering the dispersion relations 
implied by (1.7), (4.40), one can see that the two systems are formally related by a simple 
map swapping L 1/T and single-particle energies with momenta 10 : 


(k(z), e(z)) —> {-ie{z),ik(z)) 

L —>■ 1/T 


( 8 . 1 ) 


10 Notice that the identification of ik(z) as the free energy carried by a Bethe root of rapidity z is motivated 
by the discussion of Section 4.5. 
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Applying (8.1) to (5.5)-(5.7), we find a simple set of functional relations: 


p+(^)p^(2) - p+(z)p^(z) 

= 2 i sin(Lk( 2 :)), 

(8.2) 

p;(^)p v _(z)-p;(^)p v _(^) 

= —2 i sin(Lk(z)), 

(8.3) 


+ ^%)e^ x)L , a, b = +. ~, 

(8.4) 


where we require that the functions p“(z) live on a two-sheeted Riemann surface, while 
the f^s are entire. Relations (8.2)-(8.4) can be seen as a set of Baxter-like equations for 
the Hubbard Hamiltonian. For L an even integer, they admit many solutions where the f’s 
are polynomials; consequently, the p’s can be written as polynomial functions of x(z) and 
x(z), fixed in terms of their zeros on two sheets. For completeness, we can also consider 
exponential prefactors of the same type as the ones in (5.11), relabeling B —> a, n —»• /?. 
We then find easily the quantisation conditions for the following Bethe parameters: the 
zeros Sj of the p’s, and A i of the f’s. Setting k(.sy) = kj , and following the same route of 
Section 5, we find 


e ikjL+i(a —$) 

- 2 i$ TT ( A - sin(fcj) - iu \ 
M \A/ -sm(kj) + iu) 



sin (kj) — in 
sin (kj) + iu 





rriy^l 


.A'772 2zu 
A m + 2iu 


ZG{ 1 ,...,M}, 


(8.5) 

( 8 . 6 ) 


which is precisely the BA diagonalising the Hubbard Hamiltonian on the L-site chain with 
twisted boundary conditions [58, 59]: 


C L+ i,t 


ci, t e i(a -^- w) , 


C L+ 1,4. 


c u 


(8.7) 


The relation between the two Bethe Ansatz systems discussed above is directly connected 
to the path integral approach to the thermodynamics. Adopting the notation of [60], the 
latter is based on rewriting the partition function, 

Z\T> quantum = try phys (8.8) 

of a spin chain of L sites with Hamiltonian T~L, as 

■^lD quantum = Jim ^2P classical(A) A, u) 1 , (8-9) 

Af^OO TN 

where Z 2 d ciassicai(Ai) L, u ) is the partition function of an appropriately defined two-dimensional 
statistical model living on a L x N lattice. Its column-to-column transfer matrix is the 
Quantum Transfer Matrix Tqtm^)- By a 7r/2 rotation on the lattice we switch from (8.8) 
to the description 

•^lD quantum fhll tr\/ Trotter (TQ r p]y[('u)) I 1 • (8.10) 

N-¥ oo tn 

The exact Bethe Ansatz (1.14),(1.15) describes the spectrum of Tqtm m the Trotter limit 
N oo. It is an interesting open problem whether one can define the QTM directly in 
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this limit, and give it a meaningful physical interpretation as a continuum model living on 
a space of size 1/T. 

A somehow similar problem has arisen in the context of AdS /CFT integrability, where 
Zamolodchikov’s ideas on the TBA for Lorentz invariant scattering theories [61] were 
adapted to the study of the non-relativistic string sigma model dual to planar A7=4 SYM 
by considering its doubly Wick rotated counterpart (the mirror model ) [31, 62], Thanks 
to the knowledge of the action for the AdS$ x S 5 string sigma model, it has recently been 
possible to identify the corresponding mirror model as a string theory living on a mirror 
background [63]. 

Finally, we would like to mention that is also possible to introduce a mirror version 
of Takahashi’s TBA, in such a way that it is equivalent to the finite-size BA (8.5),(8.6). 
Apart from minor subtleties n , the mirror TBA equations can be obtained from (2.5)-(2.7) 
and the energy formula (2.12) through the following formal map : 

e(z) — » ik.(z), k(z) —> — ie(z), T —>■ 1/L, {B, fi} -* {ia,ifi}, (8-11) 


and by modifying the integration contours as (cf [62, 28, 29, 30]) 


/■1-M0+ r-l+i0+ roo+iO+ 

/ ^ + 

' — l+i0+ J— OO-H0+ 4l+i0+ 


( 8 . 12 ) 


so that the roles of the mirror and magic sheets of Figures 3, 4 are interchanged 12 . From 
these equations, one could in principle repeat the reduction presented in this paper and 
obtain the system (8.2)-(8.4). 


9. Conclusions 

The one-dimensional fermionic Hubbard model is one of the most interesting systems of low¬ 
dimensional condensed matter physics. Since its appearance in 1963, it has been intensively 
studied by means of exact and perturbative methods, greatly advancing the understanding 
of the physics of electron transport in ID solids. A partial grasp about the huge number of 
results on this model can be obtained by consulting [65], the book [39] and the collections of 
works in [66, 67]. The purpose of this article is to add a little piece to the jigsaw, by recasting 
the Thermodynamic Bethe Ansatz equations of Takahashi as a nonlinear Riemann-Hilbert 
problem, reminiscent of the Quantum Spectral Curve formulation recently obtained for the 
study of anomalous dimensions in AdS/CFT [32], One of the main results presented in this 
paper is a new set of nonlinear integral equations describing the thermodynamics of the 
system. In their region of validity (see discussion at the end of Section 7.2), these equations 
can be integrated numerically with very high precision and, even when implemented on 
Mathematica, the iterative algorithm converges in only a few seconds of CPU time. The 
complexity of this formulation is comparable to the system of nonlinear integral equations 

11 For example, we think the —2u/T term in (2.7) should be dropped. 

12 Due to this change in kinematics, the mirror TBA described here has some features in common with the 
TBA for the B model of [64], although the latter has a different dispersion relation. It would be interesting 
to investigate the relation between the two. 
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derived by Jiittner, Kliimper and Suzuki in [16]. However, as a consequence of the fact 
that the equations proposed here are defined on a finite support, we think that they may 
prove more convenient for the study of finite temperature correlation lengths. 
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A. The magic sheet 

In this Section we will show that the Y functions have the following cut structure on the 
magic sheet: 

• Yi^ n {z) and Y n ^(z) for n > 2 have only four branch cuts: 

z G (—1,1) ± inn , z G (—1,1) ± iu(n + 2), 

• Y].i (z) and > 2 , 2 ( 2 ) have only three branch cuts: 

z G ( — 1,1), z G (—1,1) ± 2iu. 


As we will argue in Appendix B, a stronger property is true, namely the Y functions have 
no branch points outside the positions specified above, on any Riemann sheet. However, 
here we will present only the proof on the magic sheet, which is easier as it relies only 
on the structure of the Y-system and discontinuity relations. Before we start, we need to 
rewrite the Y-system (3.4) - originally defined on the mirror section - on the magic sheet. 
Using (3.1)-(3.3), it is simple to prove that the magic-sheet version of the Y-system is 


W,t 1] W,n 1] = (1 + n,n+l)(l + Tl.n-i), n > 2 

^[ +1 ]—-pf]- = (1 + l/>n+i,i)(l + 1/Y n _gi), n > 2, 

Y n, 1 Y n, 1 

rft 11 = / i + yj, 2 \ 

y 2 [ - 1] U + l/W' 


(A.l) 

(A.2) 

(A.3) 


We will now test the cut structure of the Yi,n+i functions for n > 1. First, let us introduce 


X n = 


1 + Y i‘ 


n+ 1 ]' 


1 j_ yl n 1 

1 + 1 l,n+l 


Y 


'n+l] 


n > 2. 


1,71 


Y 


(A.4) 


1 ,n+1 
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These combinations are useful since, applying (A.l) in two elementary steps, it is possible 
to prove that 


X 


[n+3] 


1,77 


[n+2] u _ 5 


= X n X 


(i+nr-i) 


n+1 


1 + Y l 


n+3 


X 


[n-l] 


n > 2 . 


1,77+2 


(A.5) 


Notice that the term in brackets in (A.5) has no cut because it falls into the analyticity 
strip. Therefore, the cut structure of the lhs of (A.5) depends on the X n factors. We shall 
now show that X n have no cut on the real axis, for all n > 2. Starting from n = 2, we can 
compute 


v _ (1 + W 3 J)(1 + Yi$) _ Y$( 1 + Y$) _ (Y$ + 1 + Y$) _ 

^*-2 roi roi roi roi Y ' 


d2]> 


dX 


d3]> 


dl] 


dX 


Y\ ,3 


y[3]y[2] 
1 1,2 1 1,3 


Y 


[ 2 ] 


1,3 


X 


[ 2 ] 


1,3 


(i+W'i) 


. (A. 6 ) 


This expression manifestly does not have branch points on the real axis, as all terms on 
the rhs fall into their respective analyticity strips. A very similar calculation shows that 


X, = 


(i + y$)(i + xfj) _ y$(i + y/ 4 j) _ y$ + (1 + W2)(i + W 3 J) 


vWy[3] 
1 1,3 1 1,4 


= 1 + 


yP]fy[2] _i_ -I _i_ y[3]\ 
1 1,4 A1,3 “T 1 ~r J l,2l 


W 3 J(1 + ^) 


[3] N 


- [2] i 


r[3h 


= IT 


Y 


[i] 


1,4 


Y 


[i] 


1,4 


(A.7) 


(i + yPj)(i + xfjxi + yg) (i + yg) ( i + yPj)’ 

and, literally by repeating this calculation with shifted indices, one finds the general case: 


X n +1 — 1 + 


X 


[n—2] 


1,71+1 


1 _i_ y[ n 1] 
A "T 1 1,77+2 . 


1 

X 


X 


[n-l] 


1,71+2 


,1 + X/’ 


n > 2 . 


(A. 8 ) 


77+3 / 


Taking into account the analyticity strips of the functions on the rhs, this equation shows 
by induction that all X n are free of cuts. Therefore, the lhs of (A.5) is analytic in a 
neighbourhood of the real axis, meaning that X jn with n > 2 has possibly two cuts in the 
upper half plane at Iiri(z) = u(n ± 1), but no cut at Im(z) = u(n + 3). Moreover, because 
of the Y-system, no further cuts are possible and the X+ n ’s with n > 2 have only two short 
branch cuts in the upper half plane. By symmetry, the argument can be repeated for the 
lower half plane and for the Xi,i functions with n > 2 . 

The argument presented above (contrary to the remaining part of this Appendix) can be 
straightforwardly adapted to the horizontal wing of the AdS/CFT Y-systems. 

Let us now prove that X.i and Y 22 have only three cuts. Using (A.l), we can compute 


W _ v [2] 


Y£l = Y, 


2,2 


' 1 + X [ ;j ' 

,i + Vx 2 [ l 


_ y[2] 

— ■* 2,2 


' l + q + ygyy^ 
> i + (i + i/yJ3)yW < 




x 


[ 2 ] 

2,2 


yWyW 

1 1,2 ^ 2,1 


' i+n [ ,2 ' 
,1 + 1 /xj\ ] . 


R = 


y[2]y[2] 

^ 1,1 ^ 2,2 

y[l]y[l]y 

y i,2 y 2,i y 2,2 


R, 


(A.9) 

(A.10) 
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where R manifestly has no cut on the real axis and is defined by 


( l + ifj/a+yg) \ /1 + a + ifjvus \ 

“ U (dlt-vu + i/bb/ “ li + d + i/OnJ' ( ’ 

Because Yi t i(z)Y 2 1 2 (z) = e _4< ^ has no cuts outside the real axis, the only discontinuity 
can come from the factor: 


H = 1/(Y^}y 2 , 2 ). (A.12) 

Using the discontinuity relations (3.1)-(3.3), we find 

h _ i +1 Am i + i/y 2 ,2 y 2 , 2 _ , 

H 1 + Y 2 ,2 1 + Yi,i y 2 2 ' ( ' ' “ 

This shows that Y\ \ (and therefore also Y 2 2 ) does not have a branch cut with Im(z) = 4u. 
The Y-system, together with the results already obtained for the other Y functions, imply 
that no further cuts are possible. 


B. Monodromy properties of the P functions 

The purpose of this Appendix is to provide a proof for some statements made in Section 4. 

1) First, we prove that it is possible to choose the T gauges in such a way that equation 
(4.7) holds, namely 

T; i0 (z) = e~ 2 ^ z \ T^(z) = e 2 ^ z) . (B.l) 

In the proof we will use the resolvent parametrisation of Section 4.3, which, rigorously 
speaking, is valid only for the ground state. However, we expect that the result holds 
in general. 

2) Secondly, we derive the constraints (4.14), (4.15). 

3) Finally, we discuss how to infer that the second-sheet evaluation of the P functions, P, 
do not have other branch cuts apart from z£ (-1,1). This shows that the P’s live on 
a Riemann surface with only two sheets. 

To prove these properties, we will adapt many of the arguments of [44] to the present case. 


Proof of equation (4.7) for T^ 0 and Tq 1 

To verify the statement 1), let us start from the gauge T H (we restrict to the horizontal 
wing since the situation is clearly analogous for the vertical wing). This gauge is defined 
by the parametrisation 

T" n = sinh (£(! + £h + ” ] - ^h"” 1 )) /sinh (B), n> 1, 

T 0 H S = 1, seN, T 2 H i = T?J +l] l > 2. (B.2) 
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with the resolvent density p H fixed uniquely, in terms of the Y functions, through equation 
(7.1). We shall prove that 13 


TUz) = 


sinh(-Bp H (z)) 

2 sinh(2c/>(z)) sinh(i?) 


=,-2<HU 


(B.3) 


We start by rewriting some of the discontinuity relations (3.1)-(3.4) in terms of the T H 
functions. Using the properties of this gauge, equation Y\ \ = 1/Y2,2 can be written as 


( < N 
/q~H I 2,1 

'1,0 — '3,2/ '2,3 I yfT 

\ 2,1 , 


(B.4) 


while the relation Y\p 12,2 = e 4< ^ becomes 


g~ 4 0 _ 


q~H / q~H 

'1,0'2,3 _ '1,0 / '2,1 


(B.5) 


The equation involving Y \ t 2 is automatically satisfied by the parametrisation (B.2), while 
the equation for Y 2 ,i will be used later. Now, we compare the T-system equations at the 
nodes (1,1) and (2,2). On the magic sheet, these two equations read 


/■pH [H - 1] /-pH [—1] _ /yH '■pH I /-pH 

' 1,1 ' 1,1 — ' 1 , 0 ' 1,2 “T / 2,1 ? 

'"p H [+1] /-pH [—1] r-r- H [+1] /-pH [— l] /-pH ''pH '"pH I '-pH '-pH 

'2,2 '2,2 — '1,1 '1,1 '2,3 — '2,3'2,1 “T '3,2 7 1,2* 


(B.6) 

(B.7) 


The identity (B.7) was derived using T£ s = ^ . Eliminating 7j[ from these 

equations shows that 

disc 


Ti H i [+1] r ^- 11 = 7^ (r^o - 773,2/^3) = 7^ 2 r^o (1 - e 4 ^) 


(B.8) 


where we have used (B.5) in the last step. Using (B.2) in (B.8), we find the result quoted 
in (B.3). To construct the gauge T H , we can now introduce the gauge transformation factor 
h H through equation (4.25), so that h H h H = sinh(2</>sinh(£>))/sinh(7>p H )- Defining the T H 
gauge through 

= h[f ] h[r s] 71% SGN+, Tg n = 1, n G N, (B.9) 

we finally arrive at 

T^ 0 (z) = h H (z)h H (z) T? fi (z) = e~W z \ (B.10) 

13 The analogous expression for the upper wing is 

T v / \ sinh (fjpv(z)) 2 </>( z ) 

0:1 2 sinh(20(«)) sinh(/t) 
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Proof of the constraints (4.14),(4.15) 

Statement 2) can be proved by revisiting the derivation of (B.8) given above in the case of 
the gauge T H . This leads to 


= T? >2 T? i0 (l-e 4 ^). (B.ll) 

The constraint (4.14) follows from (B.ll) using (B.l) and the expression of the T H functions 
in terms of P H ’s. One can also derive an equation analogous to (B.ll) with H —> V, (/> —> — </>, 
and prove (4.15) by the same method. 


disc 


rjpH [+1] ripH [ 1] 

1,1 1,1 


Proof that the P functions have only a single cut on the second sheet 

Statement 3): let us now sketch the proof that the P’s have no cut outside the real axis 
even on their second Riemann sheet, and, therefore, live on a two-sheeted surface. The 
proof is very similar to the ones used in a more complicated context in [33, 37]. To start, 
we observe that, due to (B.3), the ratio in brackets in (B.l) and (B.5) is actually one, so 
that T 21 = 7JV Since T 2 1 = h^ ^ h H h H we find that also x has no cut on 
the real axis: 

T^ = f^. (B.12) 


Furthermore, notice that we have not used yet the discontinuity equation (3.3) involving 
Y 2 ,i. Rewriting this relation in terms of T functions in the gauge T H , we find 


rjriH [+1] rrpH [ 1] 

u 1,1 -M,l 


no T 3,i 


H [+1] Tf H [+1] 
^ 2,0 


T 


H [+ 1 ] 'jp H 1 1 


2,2 


2,2 


^3,2 T H [+1] rjpH [+1] 
3,1 1 2,0 


(B.13) 


The ratio on the lhs side can be simplified using T 2 2 = 2 \ and leads to 


1 = 


nrH upH r TP H [+1] nr H [+1] 

11 1,0 ^ 2,3 3,1 ^ 2,0 


T 3 H ,2 


= e" 4 ^ 


rrplH [+1] rjpH [+1] 

11 3,1 ^ 2,0 


rjp H [+1] rjpH [+1] 

3,1 ^ 2,0 

rjpH [+1] rrpH [+1] 

^ 3,1 ^ 2,0 


(B.ll) 


Using the T-system and (B.l), it is simple to compute the factor involving T 2 0 , and finally 
(B.14) reduces to the condition 


rjp H [+1] _ rjpH [+1] 

11 3,1 — 11 3,1 


(B.15) 


Furthermore, using the discontinuity relations 

= (i+i/y'r 1] ) / (*+v^r 1] ) > n G N+ ’ ( B - 16 ) 

which follow from the Y-system, it is possible to generalise (B.12) and (B.15) to 


disc 


rrpH \n 2] 

l n,l 


= o, 


n 


2,3, 


(B.17) 
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The set of equations (B.17) can be used to prove that the P’s live on a two-sheeted Riemann 
surface. By definition, these functions have a single cut on their first sheet, and parametrise 
some of the T H ’s as 


= p^ +s lp^[ dpHj+d, 8 € N +_ (B.18) 

Hirota equation (4.2) allows to extend (B.18) and compute any T H function in terms of the 
P H, s. However, since (4.2) is defined on the mirror section, the computation may require 
to explore the values of the P functions on the second sheet. In particular, it can be 
easily checked that T^ n with n > 1 involves (linearly) the values of P([)“" 2 \ so that 
imposing the constraints (B.17) gives information on the possible cuts of the P H ’s. For 
instance, from equation (B.12), one can recover the constraint (4.14). With the aid of this 
result, the other equations in (B.17) can be used to prove that 

pR+n] _ pH[+n] = 0j n ^ N +^ ( B . 19 ) 

Obviously, by symmetry all these results are valid both in the upper and in the lower half 
complex plane, and show that P^_ have no cuts outside the real axis. 


C. Fixing the gauge factor T H /T V 

The main purpose of this Appendix is to derive the form of the gauge transformation (4.8)- 

(4.10) relating the T H and T v functions. Along the way, we will also obtain the two NLIEs 
(C.4) and (C.8), adopted for the numerical solution method described in Section 7. 

Let us briefly discuss the logic of this Section. In Appendix B, we have encoded the 
discontinuity equations (3.1)-(3.4) in the simple statement that the P functions live on 
a two-sheeted Riemann surface and satisfy the relations (4.14),(4.15). However, some 
information is still missing to close the system. In fact, notice that, when describing the 
system with the functions P^_, there is nothing to guarantee the correct large z behaviour 
(3.5) for the Wp’s with n > 1. To fill this gap we will extract an extra equation from the 
TBA, equation (C.6), and show that it fixes the form of the gauge transformation (4.8)- 

(4.10) . Finally, we will show that the final set of equations (4.14)-(4.16) is self-consistent 
and that equation (C.6) can be rederived from them. 


Direct proof 

The Y function sitting at the central node can be parametrised equivalently using the 
gauges T H or T v . This gives the two expressions 


1 + 1/Yu = 
1 + Yu = 


sinh(20) sinh(H) Ti \ ^ 

sinh(Rp H ) 7 i H 2 

sinh(20) sinh(/t) 7^| +1 ^ \ ^ 

sinh (/ipv) 72,1 


(C.l) 

(C.2) 
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Constructing the quantity 


rM = ( l±U^M \ = ( ! + Y iM e ^ {z) \ 

W V l + Y ltl (z) ) ^ 1 + Y ltl (z) )' 


(C.3) 


and matching the alternative expressions for r(z) obtained from (C.l) and (C.2) gives the 
first equation of the NLIEs written in Section 7: 


(z)= kyyxtx = : r‘‘w 

T?,t 1] (z) T?} +1] (z) T{t 1] (z) T?} +1] (z) 


-H [-1], 


- v [+1] , 


-v[-l]. 


(C.4) 


Finding a second constraint to close the system is slightly more involved. First, using (C.l) 
and (C.2), we find: 


0 = In Yi^(z) + 4 4>{z) — In 


sinh(Bp H (z))/ sinh(5) \ _ ]n f 


sinh(—/tpv(^))/ sinh(/i) J \T7,i( z ) 


+ In 


( [+1] q~ H [—!] 

7 1,1 7 1,1 

\ t-v[+1] - 7 -v [—1] 

\ 7 1,1 ' 1,1 


(C.5) 


We can combine (C.5) with the TBA equation (2.7), simplified using the “telescoping” 
procedure discussed in Section 4.4. If the Y\^ n and Y n \ functions are parametrised using 
the T H and T v gauges, respectively, the infinite sums in (2.7) reduce to only a few terms, 
giving 


In Yqi (z) = 2u/T - 2<j)(z) - In ( ^ ] * a x (z) + In ( ^ 

\ ' 1 , 1 / \ '2,1 


(z) 


(C.6) 


By using Cauchy’s theorem and shifting the integration contours to infinity, one can write 
the last convolution term as 


In 


' 7" H 

7 1,1 


dv 


* aYz ) = - 


+ In 


disc 


2ni(v — z) 
\X} +l] (z) tY 1 \z), 


XCm T?,r‘w' 

X} + %)YJ ~%))\ 


H [ 1] / 


In 


(C.7) 


Equation (C.6) can then be rewritten as 


hi Fi,! (z) = 2u/T — 2</>(z) + j 


dv 


_i 2i\i{y — z) 


disc 


_ lu ( X l+ ‘ 1 M 7 u- 1| M ) , ]n / T "Y) 


-H [-1], 


Xl +1] (z)Y,t 1] ( z ). 




In 




(C.8) 
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Comparing (C.8) and (C.5) and using (C.4), we find 


In 


( sinh 2 (S/9 H (2:))/sinh 2 (B) 
y sinh 2 (/ipv(^))/sinh 2 (/i) 


= 4u/T + 2 


disc 


In 




Tl i [+1| (0 Tir“(v) 


•v [-1] / 


dv 

2i:i(y — z) 


— disc 


In 


vr'(z) 


Vd +1| w r lt r\z), 

Notice that the last two terms can be nicely combined as a single principal value integral: 


-v [-1]. 


(C.9) 


In 


/ sinhy PH ( z ))/sinh^) \ 4u/T + 2 C d . 
y smh 2 (/}p v (-z))/smh 2 (/i) J J- 1 


disc 


In 


' 7i"i (») -7 Ti L ~ J (^) ~ 


dv 

2tt i{y — z ) 

(C.10) 


At this point there are different ways to proceed; perhaps the simplest strategy is to apply 

- 1 dv 




w 


_i niy/l — v 2 (v — w ) 


to both sides of (C.10). Because the integrand on the rhs of (C.10) is a discontinuity of 
square root type, it can be factorised as \Jl — z 2 g(z ), where g(z) is analytic in a neigh¬ 
bourhood of the interval ( — 1,1) and therefore 


u: 


dv 


-i y/1 — v 2 (v — w)J- 


/> 


dz 


- z 2 g{z)- L - = —g(w). 


z — V 


(C.ll) 


Thus, we find 


w 


y/1 - 1/w 2 


£ 


dv 


In 


= disc 


In 


f sinh 2 (Sp H (r;))/sinh 2 (5) 
7T J_i y/l - v 2 (v — w ) y sinh 2 (/qo v (£)/ sinh 2 (fi) 

T“l + '\w) 7?}- n (w) Y 

Xi I+1| M' 7 u"‘W. ’ 


(C.12) 


(C.13) 


From the definitions of h H and h v given in (4.25), one can derive 


disc 


In 


MX) 

h v (X 


yi - ix 2 

7r 


In 


dv 


-l 


sinh(/t/9v(n)/sinh(/i) \ 

sinh(Xp H (£)/sinh(i?) J \/l — v 2 (v — z) ’ 

(C.14) 


and the comparison with (C.12) gives 


disc 


In 


MX 

MX 


= -disc 
2 




«[-i] / 


In 


(C.15) 
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Together with (C.4), this equation leads to 


disc 


In 


r; 


h[±], 


which can be interpreted as the statement that the combination 

'KAz) 


+ In 


h H (z) 

h v (z) 


= ±2 <j>(z), 


In 


= In 


niW 


+ In 


h [ + 1] ( Z )h [ H 1] (Z) \ _ [_i] [+ i] 


- (ft ^ (z) + (j)^ (z), 


(C.16) 


(C.17) 

(C.18) 


has no cuts in the whole complex plane. Using the asymptotics specified by (4.27),(4.33),(4.17), 
we conclude that the expression in (C.18) is precisely zero, proving (4.8) for the node (1,1). 
The gauge matching at the other nodes follows by the structure of the T-system. 


Converse proof 

Let us briefly discuss how equation (C.8) can be rederived from the statement that 


In 



$ ^{z) + 0[ +1 J(z) = 0, 


(C.19) 


together with the resolvent parametrisation of Section 4.3. While this may seem obvious, 
we will spell out the details in order to show that the system of nonlinear integral equations 
integrated in Section 7 are a consequence of (C.19). Just following backwards the steps in 
the derivation above and using the asymptotics liin^oo In (h H (z)/h v (z)) = —u/T, one can 
re-obtain equations (C.10)-(C.16), in particular one recovers 


/ sinh 2 (iJ/9 H (z))/sinh 2 (iJ) 
y sinh 2 (/xpv(-z))/sinh 2 (/i) 


4u/T + 2 j- disc 


In 


f 7£i 1 + 11 M7fl-‘ 1 w y 

\T^r l (v)T^t ll (v)J_ 


dv 

2ni(v — z) 

(C.20) 


Moreover, adding the two equations (C.16) shows that the two terms in (C.4) are equal. 
We can now reverse the logic and define a function Y\ \ satisfying 


f 1 + Yi,i(z) e 4< ^\ 

^ l + Yi,i(z) ) 


T"j +l \z) gCL) 

7£i HI W T“! + "(z) 


= e 4 <t>(z) 




(C.21) 


As a consequence of (C.21), Ypj automatically verifies equations (C.l), (C.2) and (C.5). 
Together with (C.20), this proves (C.8). 
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D. Dictionary 


In this Appendix, we provide a dictionary between the notation of this paper and the ones 
of [16]. First, the conventions for the coupling constant and magnetic field are related as 

u here = U rei /4, U here = ]{ rei I 10 ] /2 (D.l) 


The Quantum Transfer Matrix of [16] is characterised by three parameters: the Trotter 
number N S N, the inhonrogeneity u and the rapidity v. The Gibbs free energy at tem¬ 
perature T can be obtained by taking the Trotter limit 

/ = -T Jim \(u= -L, i/ = oV (D.2) 

iV-*oo V J 


where A is the largest eigenvalue of the QTM. The map between v and the spectral pa¬ 
rameter used in this paper, s, is 14 


1 


s = i- 


( tan 2 (u) + l ) \Z 4u2taJl2 ( u ) + ( 1 + tail2 ( w )) 2 - (D ' 3) 


2tan(i>) \tan 2 (u) + 1 / 

In particular, v = 0 corresponds to s —> oo. 

The QTM eigenvalue is given by formula (15) in [16] containing the parameters x, w, z±(x), 
z±(w), h(x), h(w) defined by 


2w = ln(tan(w)), 2x = ln(tan(u)), (D.4) 

2 h(y) = —arcsinh f —1 , z±(y) = e 2h{y)±2y . (D.5) 

Vcosh(2 y)J 

A simple calculation shows that 

z±(x) = (tan^)) 1 * 11 (1 + tan 2 (u))/ ^2u tan(u) + \J{1 + tan 2 (u)) 2 + 4u 2 tan 2 (u)^ , (D.6) 
and, using (D.3), 


z-{x) = x(s + iu), z + (x) = l/x(s — iu), (D.7) 

z^{x) — l/z^(x) = —2u±2 is, (D.8) 

Zzf(x) + l/zzf(x) = — 2(f>(s±iu). (D.9) 

The Zhukovsky map x(z) is defined in (2.13), has a single short cut and satisfies x(s) = 
— l/x(.s), and x(s) ~ 2 is at large s. Using this dictionary, equation (15) in [16] can be 
rewritten as 

'z+(x)^ 


A(u,v)/A 2 


Z-(x) 


(D.10) 


= e y+B 


+ e 


ji—B 


N/2 

Ai/A 2 + e y ~ B n 

a=l 
N/2 

A 4 /A 2 + e~ y+B n 
a=l 


s — w a — 2iu N 


Wry — s 


s — w a + 2iu N 


Wry — S 


N /T , ,, r 1 11 / s 

n l + 1/(X^- ‘(S)Zj 

3 =1 


-1 + xl 1 ] (s)/ Zj 


j-j- f 1 + l/(x[ ^(s)^') 
3 =1 


— 1 + xt +1 ] (s)/ 


14 Note that, in this Section, we use the letter s for the spectral parameter, rather than z as in the rest 
of the paper. This is to avoid confusion with the functions 2 ±(A) defined in (D.6) below. 
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The quantities A\, A 2 and A 4 are defined in [16] as 


Ai/A 2 = 
A 4 /A 2 = 


(1 - z-(w)z+(x))( 1 - z+(w)z + (x))\ N/2 


(1 + Z-(w)z + (x))( 1 + z+(w)z + (x)) t 

(1 + z-(w)/z-(x))( 1 + z + (w)/z-{x)) \ ^ 2 
(1 - z-(w)/z-(x))(- + z + (w)/z-(x))J 

A 2 = (— cos 2 (v ) cos 2 (u — u) cos 2 (v + It)) 2 ( —^ 7 — 

\z-{x)j 


N/2 




1 


z+(x)z-(w) 


N/2 


(D.ll) 

(D.12) 

(D.13) 

(D.14) 


Using the asymptotics z±(w) ~ (tan(u)) ±i = (tan(^)) ±i for N —>• 00 , it is possible to 
check that 


Jim (-1 ) N/2 A 1 /A 2 = e -Mx)+l/z+{x))/T = e 20 [ 11 ( i ) j 


N —>00 


Jim (-1 )*/ 2 A 2 /A 4 = e U-(U + 1 /*-(*))/r = e - 2 ^+H(«) } 


N—¥oo 


(D.15) 

(D.16) 


and 


lim ( cos 2 (u)e 2h ( w ' ) (1 - TzM^) ^ + 
TV —>00 V V Z-(lV) / V 


z_)_(.T)2;_(rc) 


N/2 


= e -{z-{x)-l/z + {x))/{ , n')-u/T _ p \ (0 [+1, (s)-0 [ 1] (s)) 


= ez' 


(D.17) 


Therefore, we can rewrite (D.10) in the limit of large, even N as 

A(u, u) 1 ] 00) (cos 2 (u) cos 2 (u — u) cos 2 (u T it)) 

■Q l + l/(x [+ 1 ] ( s K 0 

l=i 


N/2 


~ e 


A+s e 2 ^ [_1] + eA-^ J] 


a=l 


s — w n — 2 iu 


S — Wc 


l-xl 1 ] (s)/ Zj 


N/2 


+ e 


fi—B 


5 —2*1+1] +e -A+BJJ 




s — wTq, + 2 fu 


S — W r 


j-j- ^ 1 + l/(xl 1 ] (s)^) 
1=1 


1 — xl +1 ](s)/zj 


.(D.18) 


The factor (cos 2 (u) cos 2 (u — u) cos 2 (u + tt)) is simply a normalization (divergent in 
the N —>• 00 limit for r / 0, and converging to one for v = 0). We expect that, once this 
term is factored out, the limit N —>• 00 can be taken in the rhs of (D.18) without first 
setting v = 0 , giving a well-defined function of v (or, equivalently, of s). 

Considering the identifications discussed in Section 5, 


N/2 


F 1 (s) = e-i( B -AWu Urn TT ( 

tv-s-oo - L - L 


S - Wn 


a =1 
TV 


= e Tt^/u lim u H n (i+ i/(x(s)zj)) , 

TV-»oo 


1=1 


U = x(sj), 


(D.19) 

(D.20) 

(D.21) 
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we find a nice agreement between (D.18) and (5.35), leading to 


lim A (u,v) (cos 2 (v) cos 2 (v — u) cos 2 (v + u)) N ^ 2 = e^ + 2^ [+11 ^ 11 ( s )) Tj 1 (s). 

AT—>-oo 

(D. 22 ) 
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